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Abstract. 

We suggest a pseudospectral method for solving the three-dimensional time- 
dependent Gross-Pitaevskii (GP) equation and use it to study the resonance dynamics 
of a trapped Bose-Einstein condensate induced by a periodic variation in the atomic 
scattering length. When the frequency of oscillation of the scattering length is an 
even multiple of one of the trapping frequencies along the x, y, or z direction, the 
corresponding size of the condensate executes resonant oscillation. Using the concept 
of the differentiation matrix, the partial-differential GP equation is reduced to a set of 
coupled ordinary differential equations which is solved by a fourth-order adaptive step- 
size control Runge-Kutta method. The pseudospectral method is contrasted with the 
finite-difference method for the same problem, where the time evolution is performed 
by the Crank-Nicholson algorithm. The latter method is illustrated to be more suitable 
for a three-dimensional standing-wave optical-lattice trapping potential. 
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1. Introduction 

The experimental realization £Q of Bose-Einstein condensates (BECs) in dilute weakly- 
interacting trapped bosonic atoms at ultra-low temperature initiated intense theoretical 
effort to describe the properties of the condensate |2l El El El El IZl IHI El HOI [TH H21 H31 - 
The properties of a condensate at zero temperature are usually described by the time- 
dependent, nonlinear, mean-field Gross-Pitaevskii (GP) equation [J3]. The effect of the 
interatomic interaction leads to a nonlinear term in the GP equation which complicates 
the solution procedure. Also, to simulate the proper experimental situation one should 
be prepared to deal with an anisotropic trap [To] . 

A numerical study of the time-dependent GP equation is of interest, as this 
can provide solution to many stationary and time-evolution problems. The time- 
independent GP equation yields only the solution of stationary problems. As our 
principal interest is in time evolution problems, we shall only consider the time- 
dependent GP equation in this paper. There are many numerical methods for the 
solution of the GP equation (SHHEJEHZHSHOIIIIIIIIHIISIIIS! - 



Bose-Einstein condensation dynamics 



2 



Here we suggest a pseudospectral time-iteration method [THIE] for the solution of 
the three-dimensional GP equation with an anisotropic harmonic trap and contrast it 
with the finite-difference method [T%l ll9|. In the pseudospectral method the unknown 
wave function is expanded in terms of a set of N orthogonal polynomials. When this 
expansion is substituted into the GP equation, the (space) differential operators operate 
on a set of known polynomials and generate a differentiation matrix operating on the 
unknown coefficients. Consequently, the time-dependent partial-differential nonlinear 
GP equation in space and time variables is reduced to a set of iV coupled ordinary 
differential equations (ODEs) in time which is solved by a fourth-order adaptive step-size 
control Runge-Kutta method [20] using successive time iteration. In the pseudospectral 
method we use the Hermite polynomials to expand the wave function. In references jTT] 
pseudospectral methods have been employed for the solution of the GP equation, where a 
variable step forth order Runge-Kutta time propagator was used, as in the present work. 
In [TT] a pseudospectral fourier-sine basis was used for finite traps, and a corresponding 
complex pseudospectral basis was used for systems with periodic boundary conditions. 
However, the present study seems to be the first systematic study in dealing with the 
GP equation in three space dimensions using the pseudospectral and finite-difference 
approaches. 

In the finite-difference method the time iteration is implemented by the split- 
step Crank-Nicholson scheme. The first approach will be termed pseudospectral- 
Runge-Kutta (PSRK) method and the second finite-difference-Crank- Nicholson (FDCN) 
method in the following. Here pseudospectral and finite-difference refer to the space 
part and Runge-Kutta and Crank-Nicholson refer to the time part. However, it should 
be noted that the temporal and the spatial parts in the GP equation can be dealt 
with in independent fashions. One can combine the pseudospectral or finite difference 
approach with a specific algorithm for time stepping which can be Runge-Kutta or 
Crank-Nicholson (with or without time-splitting) scheme, or any other. In fact it is 
quite common to implement some version of the pseudospectral method with Crank- 
Nicholson time stepping when solving the Navier-Stokes equation. 

The PSRK method is illustrated by calculating the wave function and energy 
eigenvalue of the GP equation for different nonlinearity and trap symmetries in three 
dimensions as well as by studying a resonance dynamics. We find that in both cases 
the PSRK method presented here turns out to be a practical and efficient one for the 
solution of the time- dependent GP equation. 

Resonance is an interesting feature of an oscillation under the action of an external 
periodic force manifesting in a large amplitude, when the frequency of the external force 
equals a multiple of the natural frequency of oscillation. Although, the phenomenon 
of resonance is well-understood in linear problems, in nonlinear dynamics it is far 
more nontrivial. Hence, it is worthwhile to study the dynamics of resonance using 
the nonlinear GP equation. The possibility of generating a resonance in a BEC subject 
to an oscillating trapping potential has been explored previously [21]- Here, using the 
PSRK approach we study the resonance dynamics of a three-dimensional BEC subject 
to a periodic variation in the scattering length. One can have resonant oscillation in the 
x, y, and z directions, when the frequency of oscillation of scattering length equals an 
even multiple of the trapping frequency in the respective direction. Such a variation in 
the scattering length is possible near a Feshbach resonance by manipulating an external 
magnetic field [22]. There have been recent studies on this topic in one j2H] and two 
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space variables. 

Next we consider the FDCN method for the numerical solution of the time- 
dependent GP equation jUEUHllZllHl • In this approach the time- dependent GP equation 
is first discretized in space and time using a specific rule [THIHH] and the resultant set 
of equations is solved by time iteration with an initial input solution [U|3Jinj[I|. This 
procedure leads to good result in the effective one- [I] and two-space-variable jHUHllZlE] 
cases and we extend it here to the full three-dimensional case. In the two-space- variable 
case a three-step procedure jH] is first used to separate the Hamiltonian into three 
parts before applying the Crank-Nicholson scheme in two directions. In the three-space 
variable case we consider a four-step procedure where the Hamiltonian is broken up into 
four parts. We illustrate the FDCN method in three dimensions for the solution of the 
GP equation with a three-dimensional standing-wave optical-lattice potential. 

In section 2 we describe briefly the three-dimensional, time-dependent GP equation 
with an anisotropic trap. The PSRK and FDCN methods are described in section 3. In 
section 4 we report the numerical results for the wave function and energy for different 
symmetries and nonlinearities as well as an account of our study of resonance in the 
anisotropic case due to an oscillating scattering length. In section 5 we present the 
solution of the GP equation for an optical-lattice potential using the FDCN method 
and finally, in section 6 we present a discussion of our study. 

2. Nonlinear Gross-Pitaevskii Equation 

At zero temperature, the time- dependent Bose-Einstein condensate wave function 
^(r; r) at position r and time r may be described by the following mean- field nonlinear 
GP equation 0D3] 



the strength of interatomic interaction, with a the atomic scattering length. The 
normalization condition of the wave function is J dr|\I/(r; r)| 2 = 1. 

The three-dimensional trap potential is given by V(r) = ^m(u 2 x 2 + u 2 y 2 + uj 2 z 2 ), 
where u> x = u Q , u y , and uj z are the angular frequencies in the x, y and z directions, 
respectively, and r = (x, y, z) is the radial vector. The wave function can be written as 
\l/(r;T) = ip(x, y, z, r). After a transformation of variables to dimensionless quantities 

defined by x = s/2x/l, y = \^2y/l, z = y/2z/l, t = tuj , I = \J (h/muj ) and 

y, z, t) = ip(x, y, z, t)(/ 3 / the GP equation becomes 






(2.2) 
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3. Numerical Methods 

3.1. P seudo spectral Rung e-Kutta (PSRK) method 

First we describe the PSRK method for the one-dimensional GP equation in 

some detail and then indicate the necessary changes for the three-dimensional case. 
The one-dimensional GP equation is obtained by eliminating kinetic energy (derivative) 
terms in y and z, setting k = v = and eliminating the y and z dependence of in 
(ED), e.g., 



+ — +7V|0(x,t)| 2 -%■ 



0(x,t) = O, (3.1) 



dx 2 4 ' dt 

with the normalization J^° OQ dx\(j)(x,t)\ 2 = 1. 

In this method the unknown function <p(x, t) = (j)(x) is expanded in terms of a set 
of N known interpolating orthogonal functions {fj(x)}fs Q 1 as follows jTH] 

Ar ~ 1 a(x) 

<P(x) « Pn-i(x) = ^ -r-rfjWh, (3.2) 

where {xj}^, 1 is a set of distinct interpolation nodes, (f>j = (f>(xj), a(x) is a weight 
function, and the functions {fj{x)}j=o satisfy fj(xk) = 5jk (the Kronecker delta) and 
involve orthogonal polynomials of degree (N — 1), so that <p(xk) = pN~i(xk),k = 
0,1,..., AT — 1. In this work the interpolating functions {fj(x)}^S 1 are the Hermite 
polynomials Hj(x): fj(x) = Hj(x). However, one could use other polynomials, such 
as, Chebyshev, Laguerre, and Legendre. One could also consider a Fourier (spectral) 
expansion of the wave function in terms of periodic cosine and sine functions. The 
Hermite polynomials are the eigenf unctions of the linearized GP equation and hence 
already satisfy the boundary conditions of the wave function of the GP equation 
Consequently, by choosing the Hermite polynomials in the expansion, ()3.2j) satisfies the 
proper boundary conditions by construction. 

For obtaining the numerical solution, the GP equation ()3.1|) is defined and solved 
on the set of grid points Xj. The solution at any other point is obtained by using the 
interpolation formula (|3.2jl . or any other convenient interpolation rule. The advantage 
of the expansion ()3.2j) is that when it is substituted in the GP equation, the space 
derivatives operate only on the known analytic functions a(x) and fj{x), so that one 
can define a matrix for the second-order space derivative |16j : 

d 2 

(3-3) 



r>2 

k ' j dx 2 



_a(xj)' 

and the numerical differentiation process may therefore be performed as the matrix- 
vector product: Y^=q ^kj^j- Consequently, the partial differential equation (|3.1|) is 
reduced to a set of coupled ODEs in the time variable t involving fy, j = 0, 1, (N — i). 
In this way we obtain a set of ODEs by considering the original equations on a suitable 
set of discretization points (the roots of Hermite polynomials). One could have also 
used a Galerkin procedure and projected the equations, essentially by integrating 
them against Hermite polynomials for that purpose. However, we do not explore this 
possibility in this paper. 

For solving the set of ODEs we use the adaptive step-size control based on the 
embedded Runge-Kutta formulas due to Fehlberg [2H] , which gives a definite clue about 
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how to modify the step size in order to achieve a desired accuracy in a controlled way. 
For orders M higher than four of the Runge-Kutta formula, evaluation of more than M 
functions (though never more than M+2) is required. This makes the classic fourth order 
method requiring the evaluation of four functions more economic. Fehlberg suggested 
a fourth-order and a fifth-order method each requiring the evaluation of six functions. 
The difference between the results of these two gives the error 5 in the fourth-order 
method with a step size h, where 5 scales as h 5 : 5 oc h 5 . This scaling immediately 
gives the factor by which the step size h should be reduced, so that a desired 6 can be 
obtained. The detailed fourth-order and fifth-order Runge-Kutta formulas of Fehlberg 
are given in [2D]- We use these formulas with the constants given by Cash and Karp also 
tabulated in [20] • For the present problem we find that the use of Cash-Karp constants 
in the Fehlberg formulas leads to more accurate results than the original constants due 
to Fehlberg. 

Next we specialize to the case of Hermite polynomials used to generate the solution. 
Hermite polynomials are very convenient in this case as the solutions of the linear 
GP equation (J3.1j) with M = are Gaussian-type functions. They also satisfy the correct 
boundary conditions of the wave functions. In this the roots of /fjv-i(^): 

Hjsr-\{xj) = 0,j = 0,1, N — 1. The roots can be found by diagonalizing a tridiagonal 
symmetric Jacobi matrix as described in ^Bj or otherwise. The weight functions are 
taken as a(x) = exp(— x 2 /2), such that expansion ()3.2|) becomes 



7V ~ 1 exp(-x 2 /2) 
A. 

j=0 

where fj(x) are taken as 



H' N _ 1 {xj){x Xj) 

where prime refers to derivative with respect to x. Consequently, the differentiation 
matrix can be obtained from 



D 



2 



1 d 2 



exp(-x 2 /2)H N ^(x) 



(3.6) 



fc,J a(xj)H' N _ 1 (xj) dx 2 (x — Xj) 

and calculated using an algorithm described in [IS] . 

Using the differentiation matrix ()3.6|) . the GP equation is discretized. The grid 
points are the roots of the Hermite polynomial H N _i(xj) = 0. However, the actual re- 
values employed are obtained by scaling these roots by a constant factor so that most 
of the roots fall in the region where the condensate wave function is sizable and only 
a few points are located in the region where the wave function is negligible. For the 
spherically symmetric case ujq = u y = uo z , the discretization mesh in the three directions 
are identical. For the anisotropic cases, in general, the discretization points in the three 
directions are different from each other. 

Though the passage from the one-dimensional to three-dimensional PSRK method 
is formally straightforward, it involves nontrivial computational steps. The unknown 
function <f>(x, y, z, t) = <f)(x, y, z) is expanded in terms of a set of N known interpolating 
orthogonal functions {fj{x)}fj^ as follows [THj 

4>{x,y,z) ^p N ^(x,y,z) = E E E ^4 ^ fiWfMfk^ijk, (3.7) 

S i^o U> °w a (vj) a ( z k) 
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where is a set of distinct interpolation nodes, 4>ijk = 4>( x i, x j, x k), an d the 

functions a and / are defined as in the one-dimensional case, so that <p(xi,yj, = 
Pn-i( x i, Vji z k)- The differentiation matrices along x, y and 2 directions can be defined 
via ()3.6j) as in the one- dimensional case. Consequently, the partial differential GP 
equation (|2.1|) in three space variables is transformed to a set of ODEs in time variable 
on the grid points Xi, yj and Zk, which is solved by the adaptive step-size controlled 
Runge-Kutta method. The wave function at any point is then calculated using the 
interpolation formula (j3.7j) . 



3.2. Finite- difference Crank-Nicholson (FDCN) method 

The GP equation (J2.1)) has the form of the following nonlinear Schrodinger equation 

= H<P, (3.8) 

where the Hamiltonian H contains the different linear and nonlinear terms including 
the spatial derivatives. We solve this equation by time iteration after discretization in 
space and time using the finite difference scheme (TBIUHIIE]- This procedure leads to a 
set of algebraic equations. In the present split-step method the iteration is conveniently 
performed in four steps by breaking up the full Hamiltonian into different derivative 
and nonderivative parts: H = H\ + H 2 + H% + H 4 , where 



Hi = X * + + — + M^x, y, z, t) | 2 , (3.9) 

d 2 d 2 d 2 

H 2 = - ^-r, H 3 = H 4 = (3.10) 

ox 1 ay 1 oz 2 

The time variable is discretized as t n = nA where A is the time step. The solution is 
advanced first over the time step A at t n by solving ()3.8j) with H = Hi to produce an 
intermediate solution n+1 / 4 from ra , where <p n is the discretized wave function at t n . 
This propagation is performed essentially exactly for small A via 

n+i/4 = O^Hf)^ = expHAtf^", (3.11) 

where O n( ^(Hi) denotes time evolution with Hi and the suffix 'nd' refers to non- 
derivative terms. Next we perform the time propagation corresponding to the operators 
Hi,i = 2,3,4 successively via the following semi-implicit Crank-Nicholson schemes 0: 

<r +2/4 = o CN mr +1/ * = lll^l r^ (3.12) 

<r +3/4 = o CN (H 3 w +2/i = [l^/l ^ (3 - 13) 
r +1 = o CN mr +3/i = |;^g 4 ^ - +3/4 , (3.14) 

where 0(j^(Hi) denotes time evolution with Hi and the suffix 'CN' refers to Crank- 
Nicholson. Hence the final solution at time t n+ i is obtained from 

r +1 = CN (H 4 )0 CN (H 3 )0 CN (H 2 )O nd (Hi)r. (3.15) 

The details of the Crank-Nicholson discretization scheme can be found in [H]. The 
advantages of the above split-step method are the following. First, the error involved 
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in splitting the Hamiltonian is proportional to A 2 and can be neglected for small A. A 
considerable fraction Hi of the Hamiltonian is treated fairly accurately without mixing 
with the Crank-Nicholson propagation. This method can deal with a large nonlinear 
term accurately and lead to stable and accurate converged result. 

3.3. Calculational Details 

In both the FDCN and PSRK methods the time iteration is started with the following 
normalized ground-state solution of the linear GP equation (|2.1J) with M = 0: 

1/4 

exp[-(x 2 + ny 2 + vz 2 )/A\. (3.16) 

The norm of the wave function is conserved after each iteration due to the unitarity of 
the time evolution operator. However, it is of advantage to reinforce numerically the 
proper normalization of the wave function after several (100) time iterations in order 
to improve the precision of the result. Typical time step employed in the calculation 
is A = 0.001. During the iteration the coefficient n = N a/l of the nonlinear term is 
increased from at each step by Ai = 0.001 until the final value of nonlinearity n is 
attained. This corresponds to the final solution. Then several thousand time iterations 
of the equation were performed until a stable result is obtained. 

For large nonlinearity, the Thomas-Fermi (TF) solution of the GP equation is a 
better approximation to the exact result [2] than the harmonic oscillator solution (|3.16|) . 
In that case it might be advisable to use the TF solution as the initial trial input to the 
GP equation with full nonlinearity and consider time iteration of this equation without 
changing the nonlinearity. This time iteration is to be continued until a converged 
solution is obtained. However, in all the calculations reported in this paper only ()3.16j) 
is used as trial input. The use of initial TF solution did not lead to satisfactory result. 

4. Numerical Results with the PSRK Method 

4-1. Wave Function and Energy 

The present method relies on time evolution and is suitable for both stationary and 
time-evolution problems. The stationary problems are governed by a wave function 
with trivial time dependence 0(x, y, z, t) = <p(x, y, z) exp(— ijit), where [i is a real energy 
parameter. Thus the stationary wave function <fi and the parametric energy fi (the 
chemical potential) can be extracted from the evolution of the time- dependent GP 
equation over a macroscopic interval of time [HE]. Here we present results for the 
chemical potential of several BECs in three dimensions for spherically symmetric, axially 
symmetric and anisotropic cases. These traps with different geometries for 23 Na have 
been employed in experiments as well as in a time-independent solution of the three- 
dimensional GP equation [Oj. The completely anisotropic trap employs the parameters 

= 3547T rad/s, k = ^/2,u = 2 as in the experiment of Kozuma et al. [TH] . The 
cylindrically symmetric trap parameters are Uq = 33.867T rad/s, k = v = 13.585 as 
in PSj- The spherically symmetric trap parameters are Uq = 87 rad/s, k — v — 1 |27j . 
We employ the scattering length of a = 52ao of Na, where ao is the Bohr radius 

In all calculations reported in this section we used 21 Hermite polynomials each in 
x, y and z directions so that we shall be dealing with a wave function in the form of 



[x,y,z) 



K1S 



Bose-Einstein condensation dynamics 



8 



Table 1. The chemical potential /i corresponding to spherical, cylindrical and 
anisotropical geometries, respectively, for various numbers of condensate atoms iVo = 
2 q : present (f), reference {%)■ 





Spherical 


Cylindrical 


Anisotropical 


q 


n 


Mt 


n% 


n 


A*t 


n\ 


n 


Mt 


nX 


10 


0.50 


1.82 


1.825 


0.55 


17.39 


17.384 


1.77 


3.55 


3.572 


n 


1.00 


2.05 


2.065 


1.10 


19.32 


19.392 


3.54 


4.32 


4.345 


12 


1.98 


2.42 


2.435 


2.19 


22.27 


22.359 


7.09 


5.39 


5.425 


13 


3.97 


2.95 


2.970 


4.38 


26.66 


26.620 


14.18 


6.86 


6.904 


14 


7.93 


3.69 


3.719 


8.77 


32.41 


32.682 


28.35 


8.83 


8.900 


15 


15.86 


4.71 


4.743 


17.54 


41.61 


41.055 


56.70 


11.55 


11.572 



a cubic array of dimension 21 x 21 x 21. The maxima of x, y and z in discretization 
were chosen consistent with the trap parameters. Typical maxima |x|max, Mmax, and 
|^| max are °f the order of 8 for the spherical and anisotropical cases, although a smaller 
Mmax and Mmax (~ 3) together with a larger |x|max (~ 10) have been used in 
the axially symmetric case because of large distortion of trap parameters in that case 
(k = v = 13.585). 

In table 1 we list the chemical potentials for different symmetries obtained from 
the PSRK approach as a function of number of atoms iVo = 2 q in the condensate and 
compare them with those of a calculation based on a discrete variable representation of 
the time-independent GP equation 9J. In the present calculation the chemical potentials 
are extracted from results of time evolution of the GP equation. The results were 
evaluated at a space point near the center of the BEC and averaged over several samples 
of calculation. The error (standard deviation) of the time averaged chemical potential 
is typically of the order of 0.2%. Considering that the present approach is based on 
time evolution, the precision is quite satisfactory — about less than a percentage point 
of discrepancy when compared to results of 9J. It was most difficult to obtain good 
convergence in the cylindrical case with highly distorted trap (k — v — 13.585). A more 
carefully chosen values of the maxima |x|max (~ 10), and Mmax = Mmax (~ 3) were 
needed in this case. 

In the following we study the convergence of the PSRK method in the anisotropical 
case with q = 12 calculated with |x|max = 9, Mmax = 6.5 and |z|max = 5 which 
seems to be optimal and were found after some experimentation. These values were 
taken to be roughly five times the root mean square (rms) sizes (x, y, z)rms in the 
respective directions. Here we consider the convergence of these rms sizes. After the 
proper nonlinearity is introduced in the GP equation, if we continue to integrate using 
the Runge-Kutta ODE solver routine the wave function and the rms sizes fluctuate 
a little. This type of fluctuation is common to time stepping methods for a partial 
differential equation. To quantify this fluctuation we calculate the rms sizes over 250 
successive samples generated after 20 time steps of 0.01 each in the Runge-Kutta ODE 
routine. This corresponds to a total time interval of 250 x 20 x 0.01 or 50 units. Then we 
calculate the mean rms sizes and the standard deviations which we show in table 2 for 
different number of expansion points N in ()3.2j) . The numerical error in the rms sizes 
is typically less than one percent for iV > 10 and the convergence is quite satisfactory. 
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Table 2. The convergence of the rms sizes (x,y, z) rm s calculated with the PSRK 
method for the anisotropical q — 12 case for various number of space discretization 
points N. 



N 




(x)rms 


(y)rms 


(^)rms 


7 


1 


.767 ± 0.029 


1.295 ± 0.028 


0.981 ± 0.015 


8 


1 


.707 ±0.025 


1.277 ±0.018 


0.959 ±0.015 


9 


1 


.654 ±0.022 


1.251 ±0.010 


0.957 ±0.004 


10 


1 


.690 ±0.015 


1.257 ±0.008 


0.959 ±0.006 


11 


1 


.701 ± 0.002 


1.263 ±0.005 


0.960 ±0.007 


12 


1 


.691 ±0.015 


1.263 ±0.005 


0.960 ±0.007 


13 


1 


.689 ± 0.004 


1.262 ±0.012 


0.961 ±0.008 


14 


1 


.691 ±0.003 


1.262 ±0.013 


0.962 ±0.008 


15 


1 


.693 ± 0.006 


1.263 ±0.013 


0.961 ±0.009 


16 


1 


.692 ±0.005 


1.263 ±0.013 


0.962 ±0.009 



Considering that we are solving a partial differential equation in four variables this 
error is small. It will be difficult to obtain similar precision in the FDCN method in this 
problem even with a significantly larger number of space discretization points (N > 100) 
in each direction. The FDCN method is no match to the PSRK method in this problem 
with a smooth potential. However, the FDCN method seems to be very suitable for 
a rapidly varying potential, such as the optical-lattice potential considered in section 
5, which requires a large number of equally distributed spatial discretization points for 
a faithful reproduction of the potential. Although, the present PSRK method yields 
satisfactory result for the stationary problem, its main advantage lies in its ability to 
tackle time-dependent problems as we shall see in the following. 

4-2. Resonance Dynamics 

The appearance of a resonance in the oscillation of a BEC due to a periodic variation 
of the scattering length has been postulated recently in spherically and axially 
symmetric traps [21]. Here we extend this investigation to the more realistic and 
complicated case of an anisotropic trap. 

In the case of a damped classical oscillator under the action of an external periodic 
force a resonance appears when the frequency of the external force is equal to or a 
multiple of the natural frequency of oscillation of the system. The natural frequency 
of oscillation of a trapped three-dimensional BEC along a certain direction (x, y, or z) 
is twice the trapping frequency in that direction [BIEBI- For a trapped BEC subject 
to a nonlinear external force due to a periodic variation in the scattering length given 
by n = nosin(fit), in analogy with the damped classical oscillator a resonance in the 
oscillation of the BEC along a particular direction is expected when Q equals the natural 
frequency of oscillation in that direction or a multiple of this frequency. In the present 
numerical study of the anisotropic case we find that this is indeed the case. 

To investigate the phenomenon of resonance we solve ()2.1|) with TV = 
8w\/2n sin(f2£) with n = 0.3. Actually, n has to be less than a critical value n CT ^ 
(~ 0.55) in order to avoid collapse for attractive interaction (21IZ1- The actual critical 
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Figure 1. The rms sizes (cc)rms, (y)rmsj an d (z)rms y S- time for a BEC in a harmonic 
trap with k = \/2, ^ = 2 subject to a sinusoidal variation of nonlinearity n — 0.3[sin(fit) 
with (a) n = 2, (b) n = 2y% and (c) = 4. 



value depends on the asymmetry parameters k and v of the harmonic trap. For the 



spherically symmetric case k 



1. 



n 



crit 



0.575 j2J. To study the resonance 



dynamics we use 31 Hermite polynomials in each of x, y, and z directions. The maximum 
values of \x\, \y\, and \z\ employed in space discretization are each 15. The resonance 
is best studied by plotting the rms sizes (x, y, z)rms vs - time. In this study we employ 
the trap parameters uj x = ujq = 3547T rad/s, k = v2 and v = 2 [3E3> so that the 
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unit of time is oJq 1 = 0.9 ms and of length is //a/2 = 1.115 /im. For values of fl off 
the resonance the rms sizes exhibit oscillation of very small amplitude. For resonance 
frequencies fl, the rms sizes execute oscillation of large amplitude. 

To illustrate the resonance we plot in figures 1 (a), (b), and (c) (x,y, z) Tms vs. 
time for Q — 2, 2^/2, and 4, respectively. These values of correspond to the 
natural frequency of oscillation of the condensate along the x, y, and z directions, 
respectively [E1I2E]- in addition, Q = 4 is also twice the natural frequency of oscillation 
of the condensate along the x direction. So for f2 = 2 and 2\/2 the rms values of x 
and y execute resonance oscillation as shown in figures 1 (a) and (b), respectively. At 
resonance in a particular direction, the corresponding dimension increases with time in 
a oscillatory fashion. The rms values of the other components do not show resonance. 
For Q — 4, rms values of both x and z exhibit resonance. In case of (z)rms this 
corresponds to the lowest harmonic and for (x)rms this corresponds to the first excited 
state. However, the behavior of (x)rms an d (^)rms are different in this case. 

5. Numerical Results with the FDCN Method 

For a smooth trapping potential the PSRK method discussed in the last section yields 
excellent result with a smaller number of spatial discretization points which are unevenly 
distributed (at the scaled roots of the Hermite polynomial) compared to the FDCN 
method employing a relatively large number of evenly distributed spatial discretization 
points. For smooth potentials the CPU time in the PSRK method could be even an order 
of magnitude smaller than that in the FDCN method. However, the FDCN method has 
advantage in the case of a rapidly varying trapping potential requiring a large number of 
evenly distributed spatial discretization points for a proper description of the trapping 
potential. One such potential is the optical-lattice trapping potential recently used in 
BEC experiments in one [2E] and three (SOj dimensions. 

The optical-lattice potential created with the standing- wave laser field of wavelength 
A is given by V opt = V E R £f =1 sm 2 (k L Xi), with E R = h 2 k 2 L /(2m), k L = 2ir/\, and V 
the dimensionless strength of the optical-lattice potential governed by the intensity of 
the laser |3U]. In terms of the dimensionless laser wave length Ao = and the 

dimensionless standing-wave energy parameter E R /{huj) = 47t 2 /Aq, V^p^ is given by 



In the actual experiment jSD] this dimensionless standing-wave optical-lattice potential 
is superposed on the harmonic trapping potential of ()2.1|) and we present the solution 
of (|2.1|) under the action of spherically symmetric harmonic (k = v = 1) as well as the 
optical-lattice potential (|5.1j) . 

To calculate the wave function we discretize the GP equation with time step 0.001 
and space step 0.1 along x, y, and z directions spanning the cubic region — 4 < x, y, z < 4. 
Consequently, the wave function is defined on a cubic array of dimension 80 x 80 x 80. 
Starting the time iteration with the known harmonic oscillator solution for nonlinearity 
n = and Vq = 0, the nonlinearity n and the optical-lattice strength Vq are slowly 
increased by equal amount in lOOOn steps of time iteration until the desired values of 
n and Vq are obtained. Then without changing any parameters the solution is iterated 
several thousand times so that a converged solution independent of the initial inputs 




(5.1) 
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Figure 2. Three-dimensional contour plot of the interior part of the BEC ground 
state wave function under the combined action of the harmonic and the optical trap 
for n = Nqcl/1 = 10, Ao = 1 and Vq — 10 on a cubic lattice of size 3x3x3 
(—1.5 < x,y,z < 1.5): (a) view along the diagonal of the cube and (b) along one 
of the axes of the cube. 

is obtained. To illustrate the present method we calculate the wave function for the 
ground state for nonlinearity n = N a/l = 10, laser wave length Ao = 1 and optical- 
lattice strength Vq = 10. Two views of the three-dimensional contour plot of the central 
part of the wave function on a cubic lattice of size 3x3x3 are shown in figures 2 
(a) and (b) as seen from two different angles. The droplets of BEC at each pit of the 
optical-lattice potential can be identified in figures 2. There are about 10 occupied cites 
in each of x, y, and z directions of which the central part is shown in figure 2. 

6. Conclusion 

In this paper we propose and implement a PSRK method ^Bj arid contrast it with 
the FDCN method [S] for the numerical solution of the time-dependent nonlinear 
GP equation under the action of a three-dimensional trap by time iteration. In the 
PSRK method the unknown wave function is expanded in a set of known polynomials 
(Hermite). Consequently, the partial differential GP equation in three space and time 
variables becomes a set of coupled ODEs in time for the unknown coefficients, which 
is solved by a fourth-order adaptive step-size control Runge-Kutta method j20j. In the 
FDCN method the full Hamiltonian is split into the derivative and nonderivative parts. 
In this fashion the time propagation with the nonderivative parts can be treated very 
accurately. The time derivative part is treated by the Crank-Nicholson scheme in three 
independent steps. Both methods lead to stable and accurate results. The final result 
remains stable for thousands of time iteration of the GP equation. 

We applied the PSRK method for the numerical study of certain stationary and 
time-evolution problems. We solved the GP equation for spherically symmetric, axially 
symmetric and anisotropic cases and calculated the chemical potential for different 
nonlinearities. The results compare well with those obtained by a time-independent 
approach j^j. The two sets of energy values agree to within a fraction of a percentage 
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point. The numerical error in the method is found to be less than one percent with 
a small number of expansion functions (N ~ 20). The PSRK method was also used 
to study the resonance dynamics of an anisotropic BEC under the action of periodic 
sinusoidal variation of the scattering length. When this period of oscillation coincides 
with the natural frequency of oscillation of the BEC along x, y, or z directions or a 
multiple thereof the corresponding rms size executes resonant oscillation j22ll2l]- Using 
the FDCN method we study the ground state of a BEC under the combined action of 
a harmonic and a periodic optical-lattice trapping potential in three space dimensions. 

The domain of applicability of the PSRK and FDCN methods seems to be 
complementary rather than overlapping. The PSRK method is more efficient and 
economic for smooth potentials where a relatively small number of expansion functions 
and unevenly distributed spatial discretization points seems to be adequate. The FDCN 
method employs a large number of evenly distributed spatial discretization points and 
is suitable for a rapidly varying potential, such as the optical-lattice potential, favoring 
such a distribution. 
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